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ABSTRACT 

Periodic outbursts are observed in many AGNs and usually explained with a 
supermassive black hole binary (SMBHB) scenario. However, multiple periods 
are observed in some AGNs and cannot be explained with it. Here we analyze 
the periodicity of the radio light curves of AO 0235+164 at multi-frequencies and 
report the discovery of six QPOs in integer ratio 1:2:3:4:5:6 of QPO frequencies, 
of which the second with period P2 = (5.46 ± 0.47) yr is the strongest. We 
fit the radio light curves and show that the initial phases of six QPOs have 
zero or 7r differences relative to each other. We suggest a harmonic relationship 
of QPOs. The centroid frequency, relative strength, harmonic relationship and 
relative initial phases of QPOs are independent of radio frequency. The harmonic 
QPOs are likely due to the quasi-periodic injection of plasma from an oscillating 
accretion disk into the jet. We estimate the supermassive black hole mass M B h — 
(4.72 ± 2.04) x 1O 8 M and the accretion rate m ~ 0.007. With the knowledge of 
accretion disk, it implies that the inner region of accretion disk of AO 0235+164 
is a radiatively inefficient accretion flow. The oscillation accretion is due to 
the p-mode oscillation of the thick disk probably excited by a SMBHB. The 
theoretical predications of fundamental oscillation frequency and the harmonics 
are well consistent with the observations. Harmonic QPOs would be absent when 
the thick disk becomes geometrically thin due to the increase of accretion rate. 
We discuss the observations of AO 0235+164 basing on the SMBHB-thick disk 
oscillation scenario. 
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Subject headings: accretion, accretion disks — galaxies: active — BL Lacertae 
objects: individual (AO 0235+164) — galaxies: interactions — radio continuum: 
galaxies — hydrodynamics 

1. Introduction 

Blazars, consisting of BL Lac objects and flat-spectrum radio quasars, are a subclass of 
AGNs whose relativistic jet is nearly aligned to the line of sight, and show extreme variabil- 
ities at from radio through optical and X-ray to high-energy 7-ray wavelengths (Antonucci 
1993; Urry & Padovani 1995; Wagner & Witzel 1995; Ulrich et al. 1997). Long term 
multi-wavelength monitorings show that the variabilities of blazars are very complex with 
large amplitude and on time-scale ranging from tens of minutes, hours and days to months 
and years at all wavelengths. The outbursts with time-scale of order of years or longer in 
some blazars are periodic (Sillanpaa et al. 1988; Liu et al. 1995a, 1997; Raiteri et al. 
2001; Qian & Tao 2004) and of particular interest as it might connect to the presence of 
supermassive black hole binary (SMBHB) at center (Sillanpaa et al. 1988; Liu & Wu 2002; 
Qian & Tao 2004; Ostorero et al. 2004), whose ultimate coalescence would generate an 
outburst of gravitational wave radiation (Thorne & Braginskii 1976) and be the main target 
of a future space gravitational wave detector, Laser Interferometer Space Antenna (LISA) 
(Merritt et al. 2005). 

SMBHBs in galactic nuclei are expected by the hierarchical galaxy formation model 
in the cold dark matter (CDM) cosmology (Kauffmann & Haehnelt 2000) in which the 
present galaxies are the products of frequent galaxy minor mergers. Active SMBHB and final 
coalescence have been suggested to be the physical origin of the peculiar radio morphologies 
in some AGNs (for a recent review, see Komossa 2003), e.g. the helical morphology of 
radio jets (Begelman et al., 1980), the interruption and recurrence of activity in double- 
double radio galaxies (Liu, Wu & Cao 2003), the X-shaped feature of winged radio sources 
(Merritt & Ekers 2002; Liu 2004), the orbital motion of radio core (e.g. Sudou et al. 
2003). Since they were discovered in BL Lac object OJ287 (Sillanpaa et al. 1988), periodic 
optical outbursts of blazars have been ascribed to the orbital motion of SMBHB (Sillanpaa 
et al. 1988; Abraham & Romero 1999; Valtaoja et al. 2000; Liu k Wu 2002; Rieger 
2004; Lobanov & Roland 2005). However, it is unclear how SMBHBs in AGNs trigger 
the observed periodic outbursts if present. The proposed scenarios in literature include (1) 
direct impact of the secondary black hole against a standard thin accretion disk (e.g. Lehto 
& Valtonen 1996)) or an inefficient accretion flow, e.g. advection dominated accretion flow 
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(ADAF; Liu & Wu 2002), (2) rotating helical jets due to the orbital motion of SMBHB 
(Katz 1997; Villata et al., 1998; Ostorero et al. 2004; Rieger 2004), (3) periodic change 
of jet orientation introduced by the Lens-Thirring precessing of a warped disk (Abraham & 
Romero 1999; Lobanov & Roland 2005). Although the detailed physical mechanisms in 
these models are different, all of them suggest a common feature of single period in the light 
curves. 

One more complication to the investigation of the periodicity in blazars is from the 
observations of low frequency X-ray quasi-periodic oscillations (QPOs) in black hole X-ray 
binaries (Hasinger et al. 1986; van der Klis 1989; Strohmayer et al. 1996) and possibly 
in Sgr A* of our Galactic center (Torok 2005), which are believed to be due to the disk 
oscillations. Black hole X-ray binaries, so called micro-quasars, are physical analogue to 
AGNs and are powered by accretion of stellar black holes of mass about 1OM . As the 
variation time-scale of accretion disk around black holes is proportional to the mass of black 
hole and the formation of jet is coupled with accretion disk (Gallo et al. 2003; Fender & 
Belloni 2004; Fender et al. 2004), the low frequency QPO with typical frequency ~ 1Hz in 
micro-quasars corresponds to a period of order of years or longer in a blazar system of black 
hole mass larger than 1O 8 M . Disk oscillations may cause a periodic change of accretion 
and the QPOs in X-ray black hole binary system (Rezzolla et al. 2003a), which could lead 
to periodic jet emission due to the jet-disk coupling. 

To understand the physical origin of periodic outbursts, it is essential to investigate 
in details the structure of periodic outbursts in both light curves and power density spec- 
trum, and the dependence of periodicity on the observational wavelengths. For example, 
the discovery of double peaks in the periodic optical outbursts of the BL Lac object OJ287 
(Sillanpaa et al. 1996) leads to the second version of the SMBHB model (Lehto & Valto- 
nen 1996), while the investigation of the relation of the optical outburst structure and the 
radio variations leads to the third version of the SMBHB scenario (Valtaoja et al. 2000; 
Liu & Wu 2002). It is generally believed that the short-term variabilities of blazars are 
non-periodic (e.g. Wagner & Witzel 1995; Ulrich et al. 1997), but the short-term periodic 
brightness flickering may have been found to superpose on the long-term periodic variations 
in some blazars, e.g. OJ287 (Wu et al., 2006) and 3C66A (Lainela et al. 1999). In order 
to give more constraints on the models for periodic outbursts in AGNs, recently we started 
a program to investigate the relationships of periodic major outbursts and minor events and 
the long- and short-term periodicities of a large sample of AGNs (Zhao, Liu, & Wu 2006). 
Here we report the results of multiple quasi-periodic oscillations (QPOs) and the harmonic 
resonant relationship 1:2:3:. . . of QPO frequencies in BL Lac object AO 0235+164. 

BL Lac object AO 0235+164 with redshift z = 0.94 is one of the most violently variable 
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object whose periodicity of outbursts in optical and radio bands has been investigated and 
reported in literature (e.g. Webb & Smith 1989; Webb et al. 2000; Roy et al. 2000; Raiteri 
et al. 2001, 2005). Roy et al. (2000) analyzed the data of the University of Michigan 
Radio Astronomy Observatory (UMRAO) at frequencies 4.8, 8.0 and 14.5 GHz from 1975 
to 1999 and discovered a period ~ 5.8 yrs at 8.0 and 14.5 GHz, which is consistent with 
the observations of a period 5.7 yrs in the optical (Raiteri et al. 2001, 2005). Ostorero et 
al. (2004) interpreted the long term periodicity of major outbursts with a helical-jet model 
and the minor events with rotation of the non-periodic inhomogeneity. However, the radio 
outbursts are double-peaked and the burst peaks do not appear as regularly as the period 
predicted (Raiteri et al. 2001, 2005), implying that the substructure of the outburst are 
complex. On the other hand, the multiple short periods though at much lower significant 
level of signal-to-noise and the poor identification of central frequency have been claimed 
in literature. Webb & Smith (1989) first identified the three peaks with three periods 
~ 2.79 yrs, ~ 1.58 yrs, and ~ 1.29 yrs in the Fourier power spectrum of the optical light 
curves. Subsequent observations marginally confirmed the periods with 2.7 and 1.2 years in 
optical (Webb et al. 2000) and with 1.8, 2.8 and 3.7 years in radio wave-bands with the 
Discrete Correlation Function (DCF) analysis method (Raiteri et al. 2001). 

Our work is based on the three radio databases of the University of Michigan Radio As- 
tronomy Observatory (UMRAO), the National Radio Astronomy Observatory (NRAO), and 
the Metsahovi Observatory. In Sec. 2, we start the investigations from analyzing the consis- 
tency of UMRAO's and NRAO's databases, by comparing the best-sampled observations at 
8.0 GHz at UMRAO and those at 8.2 GHz at NRAO and then merging the observational 
data at these two frequencies. The periodic analysis results for the combined light curves 
at 8 GHz and identification of harmonic QPOs are given in Sec. 3. In Sec. 4, we investigate 
the dependence of QPO properties on radio frequency. After estimating the mass of central 
supermassive black hole (SMBH) and the accretion rate in Sec. 5.1, we analytically discuss 
the p-mode oscillation of thick disk in Sec. 5.2. A SMBHB-disk oscillation scenario for mul- 
tiple harmonic QPOs is suggested and discussed in Sec. 6. Our discussions and conclusions 
are given in Sec. 7. We assume a flat cosmology with H = 75 Km s~ 1 Mpc~ 1 , A = 0, and 
q = 0.5 throughout the paper, leading to a source distance d ~ AAGpc. 



2. Databases and the consistency 

UMRAO has been monitoring AO 0235+164 at 4.8, 8.0, and 14.5 GHz for about three 
decades (Aller et al. 1985, 1999), while NRAO's database contains the observational data 
at 2.5 and 8.2 GHz for about 20 years (Fiedler et al. 1987; Waltman et al. 1991; Lazio 
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et al. 2001). At frequencies 22 and 37 GHz, the Metsahovi Observatory started collecting 
the flux data since 1980 (Terasranta et al. 1998, 2004). As the central frequencies of 8.0 
GHz and 8.2 GHz are almost the same and light curves at the two frequencies are best 
sampled in both databases of UMRAO and of NRAO, we investigate the consistence of the 
observations and construct a combined light curves, based on the data from both databases. 
The observations at UMRAO last about ~ 25 years with a relatively sparse data sampling 
(~ 0.1 data points per day), while those at NRAO cover a shorter time of about 15 years 
with relatively dense data sampling (~ 0.5 data points per day). To remove the effects of 
intra-day variabilities (Romero et al. 1997) and of the noises due to unevenly sampling, we 
compose two quasi-simultaneous observational light curves by averaging the observational 
data with a 10-days time-interval, which is much shorter than the variability time-scale (of 
order of months) of major outbursts but give, respectively, about one and five observations 
in each bin at 8.0 GHz and 8.2 GHz to smooth the intra-day random variations. 

The results are shown in Fig. 1. The standard error in Fig. 1 is the statistical error. A 
least square fit with a correlation coefficient r = 0.989 gives 

F 8 . = (-0.017 ± 0.014) + (1.072 ± 0.006)F 8 . 2 , (1) 

where F 8 . and F 8 . 2 are the radio flux at 8.0 GHz (UMRAO) and at 8.2 GHz (NRAO), re- 
spectively. The results show that the observations of UMRAO and of NRAO are consistent 
with each other with a linear relation but the slope is not exact unity. So we merged the ob- 
servational data at 8.0 GHz and at 8.2 GHz and obtain a combined light curve by converting 
NRAO's observations with the fitted relations Eq. (1). The combined radio light curve at 8 
GHz is given in Fig. 7, which shows that the combined light curves fully covers several major 
outbursts in details although the combined data do not cover a significant longer time. In 
the following sections, we will give the periodic analysis results and discussions, based on 
the combined light curve at 8 GHz. In order to show how efficient and to what extent the 
combination of two databases can improve the results, we also analyze the radio light curves 
at 8.0 GHz and 8.2 GHz, respectively. 



3. Multiple QPOs in the combined 8 GHz light curve 

We analyze the periodicity of the combined light curve at 8 GHz by computing the 
Fourier power spectrum with a Lomb normalized periodo-gram (Lomb 1976; Scargle 1982), 
which is designed to properly treat the unevenly sampled data. The power spectrum is 
showed in Fig. 2. There are several very prominent peaks in the power spectrum with very 
high signal-to-noise ratios. We identify a peak with a QPO if (1) it is at least 5 times higher 
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Fig. 1. — Correlations of the observations of UMRAO at 8.0 GHz (y-axis) and of NRAO at 
8.2 GHz (x-axis). The data are 10-days averaged and the solid line is a least square fit. 



-7- 



than the nearby background noise; (2) it has a Gaussian profiles; (3) it is also detected with 
other periodic analysis methods, e.g. the Jurkevich method (Jurkevich 1971) and the 
z-transform discrete correlation function (ZDCF) technique (Alexander 1997); and/or (4) 
it correlates with other QPOs, e.g. with a harmonic relationship. 

Before discussing the periodic analysis results of the combined light curve, we firstly 
investigate the effects of the combination of the two databases. We analyze, respectively, 
with the Lomb method the light curves at 8.0 GHz at UMRAO and at 8.2 GHz at NRAO, 
and give the power spectra in Fig. 3. Fig. 2 and Fig. 3 show that all the peaks in Fig. 3 
are present in Fig. 2 and the combination of two databases significantly improves the signal- 
to-noise ratio of the power spectrum, which of the combined light curve is unprecedentedly 
high. This is because the radio light curve at 8.0 GHz has relatively sparse data points 
and misses the structure information of the major bursts, while the observations at 8.2 GHz 
span a shorter observational time interval and give only broad peaks in the power spectrum 
with small coherent quality factors. The second effect one can read from Fig. 2 and Fig. 3 
is that though most of the peaks in the power spectrum of the combined radio light curve 
can be found in the power spectra of the radio light curve at 8.0 GHz and 8.2 GHz, the 
central frequency of each peak is significantly different. A prominent peak of 8.50 years 
at the periodogram of radio light curve at 8.0 GHz is absent both at the power spectra of 
the combined radio light curve and of the light curve at 8.2 GHz. Therefore, it is probably 
spurious because of the low signal-to-noise ratio due to the incomplete coverage of major 
outbursts in the sparsely-sampled light curve at 8.0 GHz and to the loss of the burst structure 
information. To solve the problem, much more intensive observations are need. 

The peaks in Fig. 2 are Gaussian, though some of them overlap with each other. There- 
fore, we fit each peak profile with a Gaussian function or, if necessary, fit several peak profiles 
together with Gaussian functions. The fitting results are given in Table 1. In Table 1, Col. 1 
gives the observational radio frequency z/ obs , the averaged period of the fundamental QPO 
and their standard errors. We will discuss the definition of fundamental QPO and how to 
calculate the average period later. The centroid frequency v and the standard error a of 
each peak are given in Col. 2 and indicated in Fig. 2. The fitting goodness x 2 is presented 
in Col. 6. The standard error a is that of a fitting Gaussian function and is considered to 
contain all the effects, including random variations in the exact quasi-periodic oscillation, 
the large and changing width of the outburst structure, the poor coverage of some outbursts, 
the random variations in intensity, and in particular the intrinsic coherence of variability. 
To quantify the coherence of the variabilities, we calculate the quality factor Q = Z//HWHM 
and give it in Col. 5, where HWHM is the half width at half maximum of a peak. The 
quality factor Q in Table 1 increases with the QPO frequency v and correlates with the 
repetitions of a QPO in the combined radio light curve, which is calculated with N = uAT, 
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Fig. 2. — The power spectrum of the combined radio light curve at 8 GHz. Fitted centroid 
periods of QPOs are indicated. The dashed and the dash-dotted lines are the la and 3a 
significance levels, respectively. 
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where AT ~ 26.3 yr is the total observation time of the combined data sample. A strong 
correlation of the quality factor Q and the repetition number N is showed in Fig. 4, implying 
that the coherence factor Q- value and the quality of QPOs are limited by the duration of the 
monitoring program. Therefore, to understand the intrinsic coherence of a QPO, we have to 
monitor the object for a much longer time. 

In the Lomb periodogram of the combined 8 GHz light curve, we identify six QPOs of 
frequencies (periods) v x = 0.2279 x 10" 3 day _1 (P l = 12.02 yr), v 2 = 0.5029 x lO^day -1 
(P 2 = 5.45 yr), i/ 3 = 0.7568 x lO^day" 1 (P 3 = 3.62 yr), v A = 0.9691 x lO^day" 1 (P 4 = 
2.83 yr), u 5 = 1.275 x lO^day" 1 (P 5 = 2.15 yr), and v % = 1.511 x lO^day" 1 (P 6 = 1.81 yr). 
The periods of the QPOs and the standard errors are given in Col. 4 of Table 1. To estimate 
the relative contribution of each QPO signal to the total flux variations, we compute the 
relative root mean square (rms) in percentage and give the results in Col. 3 of Table 1. The 
results show that the QPO of period P = 5.45 yrs is the strongest with rms ~ 10 per cent, 
though the other QPOs are also significant. We identify the frequency of the strongest QPO 
with the fundamental one with frequency u s . 

In order to confirm the detection of multiple QPOs, we now apply the least variance 
analysis introduced by Jurkevich (1971) and the z-transform discrete correlation function 
(ZDCF) technique (Alexander 1997) to the combined 8 GHz light curve. Jurkevich method 
has been successfully developed to search for periods in the optical light curves of some 
blazars (e.g. Kidger et al. 1992; Liu et al. 1995a, 1997). The method is based on the Ex- 
pected Mean Square Deviation and specifically designed for the unevenly sampled datasets. 
It tests a run of trial periods. A light curve is first folded in a trial period and the variance 
V£ is computed first for each phase bin and then over a whole period. For a false trial, V£ is 
almost constant but for a true one it reaches a minimum. Instead of the normalized itself, 
a parameter / = (1 — V£)/V£ is usually introduced to described the significance of a period 
in a light curve (Kidger et al. 1992; Liu et al. 1995a). It is suggested that / > 0.25 implies 
a strong period. A more general description and a robust test of the method was given by 
Liu et al. (1997). The analysis results with Jurkevich method are given in Fig. 5. In 
the V^-plot with very high signal-to-noise ratio, there are several very prominent minima 
with broad width at the trial periods P = 1.78 yr, 3.46 yr, 5.34 yr, 7.90 yr, 11.60 yr and 
15.10 yr. The last period may not be real, because it is very long and the data sample covers 
it less than two times. Among all the six QPOs detected with the Lomb power spectrum 
method, strong QPOs of periods Pi, P 2 , P3, and P 6 with rms larger than 7 per cent are 
clearly detected by the Jurkevich method, though the central frequencies of the periods 
shift toward to the lower end and cannot be determined accurately due to large noise. Weak 
QPOs V4 and 1/5 with rms < 6% are insignificant in the normalized V^-plot and are nearly 
drowned out by noise. In Fig. 5, a period P = 7.90 yr is significant but not present in the 
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Lomb periodogram shown in Fig. 2. One possibility is that the QPO at Pi = 12.02 yr with 
low quality factor Q in the Lomb periodogram may consist of several QPOs (e.g. those 
periods P = 7.90, 11.6 yr, and 15.10 yr in the normalized V^-plot) and cannot be resolved 
due to the very low quality factor Q. However, it is more likely that the period P = 7.90 yr 
in Fig. 5 does not represent a real QPO but just a spurious period at the position of two 
and four times of the two strong QPOs P 3 = 3.62 yr and P 6 = 1.81 yr, respectively. Or, the 
period P = 7.90 yr is just the result of incorrect identification of noise. 

We also analyze the combined radio light curve with the z-transform discrete correla- 
tion function (ZDCF) technique (Alexander 1997), which is designed properly to treat the 
unevenly sampled datasets. Discrete self-correlation function has been successfully used to 
search for periods in the light curves of AO 0235+164 (e.g. Raiteri et al. 2001) and ZDCF 
is basically the same as DCF. The self-ZDCF function of the combined radio light curve at 
8 GHz is showed in Fig. 6. In self-ZDCF, a local maximum implies a period. In Fig. 6, 
all the maxima have very broad profile and the period is difficult to determine accurately. 
Because our purpose is to confirm the periodic analysis results obtained with the Lomb 
power spectrum method and the Jurkevich method, we do not try to compute periods from 
ZDCF accurately but estimate a period at peak. The identified five QPOs have periods 
P = 11.34 yr, 7.39 yr, 5.53 yr, 3.33 yr, 1.88 yr. The periodic analysis results with the ZDCF 
and Jurkevich methods are consistent with each other within the estimate errors but not 
fully consistent with the results with the Lomb power spectrum method. Like the Jurkevich 
method, ZDCF does not detect the two weak periods P 4 = 2.83 yrs and P 5 = 2.15 yrs ob- 
tained with Lomb's power spectrum method. It might be due to the broad profiles of the 
strong periods 1.88 yrs and 3.33 yrs at ZDCF and/or to the weakness of the two periods, so 
that the two weak QPOs are absent in the ZDCF. 

The conclusions from three different periodic analysis methods are that the strong QPOs 
with periods P = 12.02, 5.45, 3.62, and 1.81 yrs are detected with unprecedentedly high 
singal-to-noise ratio in the combined 8 GHz radio light curves. Two weak QPOs with P = 
2.83 yr and 2.15 yr are detected by the Lomb power spectrum method but are not detected 
by the other two methods due to the weakness. They will be confirmed with the periodic 
investigation of radio light curves at multiple radio frequencies in Sec. 4.2. 
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4. Harmonic relationship and dependence of QPOs on radio frequencies 

4.1. Harmonic relationship of QPOs 

The ratio of centroid frequency and the fundamental frequency, f u = u/u si of six QPOs 
is given in Col. 8 in Table 1. We do not compute the error for the ratio as the standard errors 
a of the frequencies given in Col. 2 are limited by the time coverage of the data sample. The 
ratios f v of six QPOs can be given with a simple relation f v = 1/2 with I — 1, 2, 3, 4, 5, 
and 6, or all six QPOs have harmonic relation V\ : v 2 : v-i '■ : u 5 : z/ 6 = 1 : 2 : 3 : 4 : 5 : 6. 
The identification is given in Col. 9 of Table 1. With the identification, we calculate the 
averaged fundamental frequency V s and the standard error a s with P s = a s ^2{vj/crj)/N and 
l/^s = XX.? / ,(T j) /^N with j — 2, . . ., and 6, where N is the total number of QPOs used in 
the calculation. In the computation, we do not take into account the QPO with the lowest 
frequency v\ (the longest period P\ = 12.02 yr) due to the low quality factor Q = 1.76. The 
averaged frequency of the fundamental QPO is v s = (0.5014 ± 0.0434) x 10~ 3 day -1 (period 
P s = 5.46 ± 0.47 yr). With the least square method, we fit the combined 8 GHz radio light 
curve both with the fundamental QPO alone (x 2 = 0.7539) and with the six detected QPOs 
(x 2 = 0.3235), and give the results in Fig. 7. Fig. 7 shows that the fundamental QPO alone 
cannot satisfactorily reproduce the radio light curve while the fit with six QPOs can produce 
it very well until about the year 2000. The results imply that the contributions from the 
harmonic QPOs are significant, and the sub-peaks of the outbursts are not random but have 
the same physical origin as that of the major outbursts. The differences of the initial phases 
of six QPO harmonics relative to the fundamental are listed in the last column of Table 1. 
The results show that the QPOs vi, u 2 , v 3 , v 5 , and z/ 6 almost have the same initial phases 
with A0o ~ while A0o ~ 71 f° r the 4th QPOs 1S4. The initial phase difference A0o of 
the 1st and the 5th QPOs are deviated from zero because the large fit uncertainties due to 
the low quality factor Q for the 1st harmonics or to the weakness (rms = 4.70%) of the 
5th harmonics, respectively. The relative initial phases of QPOs are consistent with the 
suggestion of a harmonic relationship among QPOs. 



4.2. Frequency dependence of harmonic QPOs 

In this section, we present QPO analyses with the radio light curves at 4.8 GHz, 14.5 
GHz, 22 GHz, and 37 GHz and investigate the dependence of the harmonic QPOs on radio 
frequency. Since our purpose is to investigate the QPOs and the frequency dependence, we 
give only the analysis results with the Lomb power spectrum method. The radio data at 4.8 
GHz and 14.5 GHz from UMRAO and at 22 GHz and 37 GHz from Metsahovi Observatory 
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are plotted in Fig. 7. The periodic analysis results are given in Fig. 8. 

Fig. 8 shows that the Lomb periodograms of radio light curves at the four radio fre- 
quencies are very noisy and the signal-to-noises of the Lomb periodogram are significantly 
lower than that of the combined radio light curve at 8 GHz. Among six QPOs obtained at 8 
GHz, the five strongest QPOs are detected in the radio light curves at all the radio frequen- 
cies, while the weakest QPO of period P5 = 2.15 yr and rms = 4.70 per cent at 8 GHz is 
significantly detected only at 14.5 GHz, at which frequency the radio light curve is relatively 
well sampled. At 4.8 GHz, 22 GHz and 37 GHz, the weakest QPO is absent from the Lomb 
power spectra probably due to poor sampling and to the missing of the information on the 
minor events. Like at 8 GHz, the peak at a period with P ~ 7yr obtained both with the 
Jurkevich and the ZDCF methods, is not present in the Lomb periodograms at all the 
radio frequencies in Fig. 8. Table 1 lists the fitted centroid frequency with Gaussian function, 
fitting goodness y 2 , relative strength rms, quality factor Q, the ratio of the frequencies of 
harmonics and the fundamental QPOs, and the averaged fundamental period P s . Again, the 
QPOs clearly have a harmonic relationship and the identification of the harmonic integer is 
given in Table 1. Like what we did for the combined 8 GHz light curve, with the least square 
method we also fit the radio light curves at 4.8 GHz, 14.5 GHz, 22 GHz and 37 GHz both 
with the fundamental QPO u s alone and the six QPOs together, respectively. The fitted 
light curves are given in Fig. 7 and the differences of the initial phases of harmonic and the 
fundamental QPOs are listed in Col. 10 of Table 1. From the fitting results in Table 1, the 
differences of initial phases of the six QPOs at the four radio frequencies are consistent with 
the results obtained at 8 GHz and five QPOs vi, v 2 , v 3 , and u Q have almost the same 
initial phases but about 7r difference from that of the 4th one. The results strongly support 
the identification of six QPOs and the harmonic relationship. 

The dependence of QPO centroid frequency on the observational radio frequencies is 
given in Fig. 9. Fig. 9 shows that QPO frequencies are independent of the radio frequencies 
and the ratios of the frequencies of harmonic QPOs at the four radio wave-bands and the 
fundamental one with u s = 0.5029 x 10~ 3 day~ 1 at 8 GHz are also independent of the radio 
frequencies, implying that the QPO centroid frequencies and the harmonic relationship are 
not determined by the different radio emission regions in relativistic jet. However, the 
independence does not mean that one could detect all the harmonic components of QPOs at 
all wave-bands. To answer this question, we investigate the dependence of relative amplitude 
rms on radio frequencies. Fig. 10 shows the results that the relative amplitude rms of all 
QPOs is a function of radio frequency and the third QPO has the strongest dependence. 
However, the conclusion is very sensitive to the periodic analysis results of radio light curves 
at 37 GHz, which is poorly sampled. Therefore, the decrease of QPO amplitudes may be due 
to the loss of information on high frequency QPOs in the poorly-sampled radio light curves. 
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To eliminate the effect of the incomplete coverage of outburst structures on the analysis 
results of QPO harmonic components, we calculate the ratio of rms of the harmonics and 
rms(f s ) of the fundamental QPO. The results are given in Fig. 10. In Fig. 10, rms / rms{y s ) 
for the 1st, 4th and 6th components of the harmonic QPOs is nearly independent of frequency 
while for the 3rd harmonics the ratio decreases with radio frequency. Again, the conclusion 
is sensitive to the periodic analysis results of the 37 GHz light curve. If the analysis results 
of QPOs at 37 GHz are not considered, the harmonics of QPOs are nearly independent of 
radio frequencies. The results suggest that if one detects the fundamental QPO of period 
P s = 5.46 yr in a well-sampled light curve at high frequency, e.g. in the optical wave-bands, 
one would simultaneously detect the other QPOs with periods Pi ~ 10.9 yr, P 4 ~ 2.8 yr, and 
P 6 ~ 1.8 yr. But, the QPO of period P3 = 3.6 yr may not be observable in optical wave-band. 
Although the conclusion needs to be tested with much more observations at 37 GHz, our 
results suggest that the properties of QPOs would be independent of radio frequency and 
not be determined by the physical process in the radio emission regions of jet. 



5. Central black hole mass and thick disk oscillations 

Raiteri et al. (2001) reported the detection of a period P = 5.7 ± 0.5 yrs in the optical 
and radio light curves of AO 0235+164. Ostorero et al. (2004) interpreted it with a helical 
jet model, in which the quasi-periodic radio-optical light curve and the spectral energy 
distributions (SEDs) are caused by the orientation variation of a helical inhomogeneous, 
non-thermally emitting jet probably due to the orbital motion of the primary black hole in 
a binary system. In the helical jet model, the minor burst events are taken as non-periodic 
flux fluctuations due to the random plasma instabilities. However, our periodic analyses 
show that the major outbursts and the minor burst events are the products of combinations 
of six harmonic QPOs, which suggest a quasi-periodic radio emission from relativistic jets 
and a physical origin of harmonic QPOs independent of the radio emission regions in the 
jet. The multiplicity, the harmonic relationship, the zero differences of initial phases, and 
the independence of relative rms on radio frequency imply that the harmonic QPOs in the 
radio light curves of AO 0235+164 are most likely due to the quasi-periodic injection of 
plasma from accretion disk into the relativistic jet. Because of the jet-disk coupling (Fender 
& Belloni 2004), such a quasi-periodic plasma injection implies a quasi-periodic oscillation 
of accretion around SMBH. 

Both thin and thick accretion disks can oscillate quasi-periodically, but only the oscil- 
lation of a thick disk can be global and trigger the quasi-periodic accretion with harmonic 
frequencies. Rezzolla et al. (2003a) and Zanotti et al. (2003) found that the accretion 



-14- 



rate of a tori or thick disk of finite radial extent is quasi-periodic due to the global p-mode 
oscillation, which could be excited by a global perturbation (Rezzolla et al. 2003b; Zanotti 
et al. 2003) or a local periodic strong agent (Rubio-Herra & Lee 2005a,b). Such a torus or 
thick disk of finite radial extent could be a optically thin and geometrically thick radiatively 
inefficient accretion flow (RIAF) or advection dominated accretion flow (ADAF) between a 
central supermassive black hole and an outer truncated geometrically thin disk. It is clear 
now that such a configuration should be present in a black hole accretion system if the di- 
mensionless accretion rate m = M/M Edd is lower than a critical value m cr ~ 0.02 (e.g. Meyer 
& Meyer-Hofmeister 1994; Narayan & Yi 1995; Narayan 2002; Esin et al. 1997; Meyer 
et al. 2000; Meyer-Hofmeister et al. 2005), where M E dd = -^Edd/O.lc 2 is the Eddington 
accretion rate relative to the Eddington luminosity Ls&d — 1-26 x 10 46 erg s -1 (M/1O 8 M ). 
Because the p-mode oscillation frequency of thick accretion disk is determined by the disk 
radial extent defined with the dimensionless accretion rate, we now estimate the mass Mbh 
of central black hole and then the dimensionless accretion rate m of AO 0235+164. 



5.1. Black hole mass and accretion rate 



As there is no valid estimation of the black hole mass of AO 0235+164 in literature, 
we estimate its black hole mass and the accretion rate in this section, before we start the 
investigations of the physical origin of QPOs. Because of the high redshift (z = 0.94) and the 
very bright central nucleus, it is impossible to resolve the host galaxy of AO 0235+164 even 
with the Hubble Space Telescope and to infer the black hole mass with the tight relation 
of black hole mass and the properties of galaxy bulge (e.g. stellar velocity dispersion a). 
However, we can estimate the black hole mass with some recently suggested relations of the 
emission line properties and central black hole mass in AGNs. 

The reverberation mapping studies reveal an empirical relation of broad emission line 
region (BLR) size and the optical continuum luminosity at rest-frame 5100A (Kaspi et al. 
2005) 



R 



BLR 



(2.23 + 0.21) 



AL A (5100A) 
10 44 ergs s- 1 



0.69±0.05 



(2) 



10 It — days 

Together with the FWHM of H^ emission line, which is a measurement of the characteristic 
velocity of broad emission line clouds, one can estimate the central black hole mass of an 
AGN with the relation (Peterson et al., 2004) 



M BH = 1-936 x 1O 5 M 
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BLR 



It — days 



FWHM(H /3 ) 
10 3 km s" 1 



1 2 



(3) 
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The environment of AO 0235+164 is rather complex due to the presence of several 
foreground intervening galaxies within a few arc-seconds. Great efforts have been done to 
study the emission line features of the object with spectroscopic observations. Using the 
3m Shane telescope at Lick Observatory, Cohen et al. (1987) identified the emission lines 
of Mgll, [NeV] and [Oil] at the source red shift. With the 4m William Herschel Telescope, 
Nilsson et al. (1996) detected not only the [NeV] and [Oil] lines but also prominent hydrogen 
emission lines H<5 and H 7 at the same redshift. Nilsson et al. (1996) obtained the FWHM 
and flux (in the observer's frame) of H,5(A4102) and H 7 (A4340) emission lines to be 3600 ± 
400 km s" 1 , SOxlO-^ergs^cm- 2 and 3400±400km s~\ 97xl0~ 17 erg s^cm" 2 , respectively. 
The equivalent widths of the H<5 and H 7 lines are 2.6A and 7.2A (all the quantities are 
given in rest frame except noted). With the assumption of a power-law optical continuum 
(F u oc v~ a \ we can estimate the flux at 5100A from the flux in the bluer band. However, 
deriving an accurate value of spectral index a is not a simple task for AO 0235+164 because 
of its strong variation (Smith et al. 1987). A visual inspection of Fig. 1 of Cohen et 
al. (1987) reveals that at 4500A (observer's frame) the continuum flux density is about 
0.5xl0~ 16 erg s _1 cm _2 A _1 , while at 7500A it is about l.lxl0~ 16 erg s _1 cm _2 A _1 . From these 
flux values we can derive a=3.54, which is in fair agreement with the values of a ~ 3.6 
corresponding to a faint state of the source reported in Fig. 7 of Raiteri et al. (2001). 
However, when one takes into account the absorption by the intervening system at z=0.524, 
this value nearly halves. This can be inferred by the intrinsic (dereddened) < B — R >= 1.72 
derived by Raiteri et al. (2005) which, when considering the amount of extra absorption 
tabulated in their Table 5, gives < a >= 1.75. We notice that this a value means a relatively 
flat shape in the optical band of the flux density versus wavelength plot (Fa oc A a ~ 2 ), 
which is also consistent with the dereddenned spectrum shown in Fig. 7 of Junkkarinen 
et al. (2004). Adopting a=1.75, from the observed flux at 7500l (rest-frame 3866A) 
in Fig. 1 of Cohen et al. (1987) we calculate the flux density at rest-frame 5100A as 
7.5xl0~ 16 erg s~ 1 cm~ 2 A~ 1 , while from the continuum flux at 4102A estimated from the H<j 
line flux and equivalent width in Nilsson et al. (1996) we calculate the flux density at rest- 
frame 5100A as 5.3x 10~ 16 erg s _1 cm~ 2 A~ 1 . From the average value of these two estimates, we 
obtain the monochromatic luminosity at 5100A with ALa(5100A) = 7.48 x 10 45 erg s -1 . From 
Eq. (2), we estimate the BLR size as -Rblr = 438 ± 103 It — days. Because of no detection 
of H/3 for AO 0235+164 in literature, we try to estimate its FWHM from the observations 
of other Hydrogen Balmer lines. As the observed H 7 line is probably contaminated by 
[OIII](A4363) line, we make the estimation with H^ line data. From the observations of 50 
quasars in LBQS given by Forster et al. (2001), we derive an empirical relation between the 
FWHM values of Hg and H 5 lines with the OLS bisector method (Isobe et al. 1990) 



FWHM(H /3 ) = (415.6 ± 348.8) + (1.07 ± 0.15)FWHM(H 4 ). 



(4) 
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From the observations of the FWHM of H 5 line, we have FWHM(H /3 ) ~ 4267±772 km s" 1 for 
AO 0235+164. From Eq. (3), we estimate the central black hole mass of AO 0235+164 to be 
M BH - (1.54 + 0.67) x 1O 9 M . However, AO 0235+164 is a blazar-type object and its optical 
continuum luminosity is generally dominated by the beamed synchrotron emission from a 
relativistic jet rather by the emission from accretion disk. The southern intervening galaxy 
with z=0.524 may also partly contribute to the optical continuum luminosity (Nilsson et al. 
1996). Therefore, the estimated black hole mass with the optical continuum luminosity is 
overestimated and can be taken only as an upper limit. 

To avoid such an overestimation of black hole mass, Wu et al. (2004) argued that the 
emission line luminosity is probably a better tracer of the ionizing continuum luminosity of 
blazars and proposed a new empirical relation between the BLR size and the Hp emission 
line luminosity 



lgi?BLR(lt-days) = (1.381 + 0.080) + (0.684 ±0.106) lg(L^ /10 42 ergs s" 1 ). (5) 



Because no detection of emission line is reported for AO 0235+164 in literature, we 
estimate the luminosity from the observations of H<5 and H 7 line fluxes. Zheng (1988) 
studied the variations of Balmer decrements of several quasars. From his table 2 and table 
3, we get an averaged flux ratios 5 and 2.4 of H^ to H5 and H# to H 7 , respectively. From 
the observations of emission lines by Nilsson et al. (1996), we obtain the rest-frame 
luminosity for AO 0235+164 with Lh — 1-68 x 10 43 ergs s -1 from the observed flux of 
US emission line. The estimated luminosity from the observed H 7 flux is about 15 per 
cent higher, which might be due to the contamination of [OIII](A4363) line. From Eq. (5) 
and the H3 luminosity estimated from H 5 emission line flux, we calculate the BLR size as 
-Rblr = 166 ± 64 It — days, which is much smaller than that estimated from the optical 
continuum luminosity. With the estimated BLR size and H/3 FWHM value, from Eq. (3) we 
estimate the black hole mass of AO 0235+164 of M BH = (5.85 ± 3.09) x 1O 8 M , which is 
much smaller than that obtained with the optical continuum luminosity method. 

McLure & Dunlop (2004) suggested to estimate the black hole mass of AGNs from the 
UV continuum luminosity and the FWHM of MgII(2788A) emission line with an empirical 
relation (with cosmological parameter Hq = 70 km s" 1 Mpc -1 , Q m = 0.3 and = 0.7; the 
distance of AO 0235+164 is 6.1Gpc in this case) 



For AO 0235+164, the Mgll emission line is observed to have the FWHM of 3100 km s" 1 
and the flux of 1.24 x 10~ 15 ergs cm" 2 s -1 (Cohen et al. 1987). The equivalent width 
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and the continuum flux are 15. 7 A and 0.8 x 10~ 16 ergs cm" 2 s _1 A -1 , respectively. Using 
the observations of the Mgll emission line and the continuum, we derive the rest-frame 
continuum luminosity at 30001 of AL ° = 7.37 x 10 45 erg s -1 (using the same cosmological 

J 3000A to V to to 

parameters as in McLure & Dunlop (2004)). From Eq. (6), we obtain a black hole mass 
-^bh = 4.42 x 10 8 M Q , which is approximately the same as the estimation obtained with 
the luminosity of H/3 emission line. However, the same as the optical continuum luminosity, 
the UV continuum luminosity of A0235+164 is also contributed by the jet emission and 
probably by the intervening galaxy. On the other hand, a recent study of Junkkarinen et 
al. (2004) indicated that the UV continuum of AO 0235+16 is heavily absorbed by the 
Galaxy and the intervening materials at z=0.524, which implies that the de-reddened UV 
continuum is much larger than the observed one (see their Fig. 7). Therefore, the UV 
continuum luminosity of AO 0235+164 suffers large uncertainties and the estimated black 
hole mass derived with it may be unreliable and will be considered only as a reference. 

It has been found that there is a tight correlation between the black hole mass and the 
bulge velocity dispersion for nearby galaxies (Gebhardt et al. 2000; Merritt & Ferrarese 
2001; Tremaine et al. 2002), which can be expressed as (Tremaine et al. 2002) 

Log(M BH /M & ) = 8.13 + 0.06 + (4.02 + 0.32) Log (a/ 200 km s" 1 ). (7) 

AGNs seem to follow the same correlation as nearby galaxies (Ferrarese et al. 2001; Onken 
et al. 2004). On the other hand, there is also a tight correlation between the central velocity 
dispersion of galactic bulge and the FWHM of [OIII] (50071) emission line in AGNs with 
a w FWHM([OIII])/2.35 (Nelson & Whittle 1995). This suggests that black hole mass 
could be roughly estimated from the observed FWHM values of narrow forbidden lines in 
AGNs. For AO 0235+164, both Cohen et al. (1987) and Nilsson et al. (1996) have detected 
narrow forbidden lines at redshift of z = 0.94. Cohen et al. (1987) estimated the FWHM 
values of [Ne V] (34261) and [Oil] (37271) as 600 km s' 1 and 200A;m s~\ respectively. 
Nilsson et al. (1996) obtained the FWHM values of [Ne V] (33461), [Ne V] (34261) and 
[Oil] (37271) as 700±100A;m s~\ 400+100 km s" 1 and 600±100A;m s~\ respectively. From 
the emission line data of quasars in LBQS (Forster et al. 2001), we know that the average 
FWHM value of [OIII] (50071) is comparable with those of [Oil] (37271) and other narrow 
forbidden lines. Therefore we roughly take 600A;m s _1 as the FWHM value of [OIII] (50071). 
Using a FWHM([OI 1 1]) /2.35 : we derive the value of a central velocity dispersion a = 
255km s -1 . From Eq. (7), we estimate the black hole mass of AO 0235+164 as M BH ps 
(3.58 ± 2.65) x 1O 8 M , which is consistent with the value obtained with the H/3-luminosity. 

We have adopted four different kinds of methods to derive the SMBH mass of AO 
0235+16. Except the optical continuum luminosity method which apprently overestimates 
the SMBH mass, the other three methods give the SMBH mass in a range from 3 x 1O 8 M to 
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6 x 1O 8 M , with typical uncertainties of a factor of a few. In the following discussions, we take 
the averaged mass M BH = (4.72 ±2.04) x 10 s M for the central black hole of AO 0235+164. 
However, we must keep in mind that the value is derived either with the indirect estimated 
H/3 and [OIII](5007A) emission line properties or with the UV continuum luminosity which 
suffers from serious contaminations from the jet, southern intervening galaxy and extra- 
absorptions. In order to obtain more reliable black hole mass, further studies are absolutely 
needed to accurately determine the emission line and continuum properties of AO 0235+164. 

Wang et al. (2004) calculated the broad line luminosity Lblr — 4.1 x 10 43 ergs s _1 
for AO 0235+164 (adjusted to fit the cosmological parameters used here). If we adopted 
a typical covering factor C ~ 0.1 for the broad-emission line region (Netzer 1990), the 
bolometric or the ionization luminosity is L bol = L ion ~ L BLR /C ~ 4.1 x 10 44 ergs s -1 . 
From the estimated bolometric luminosity and black hole mass, we obtain the dimensionless 
accretion rate m~7x 10 _3 (?7/0.1) _1 of AO 0235+164, where rj is the conversion efficiency 
of matter to energy. 

5.2. Thick accretion disk p-mode oscillations and physical origin of multiple 

QPOs 

The estimated dimensionless accretion rate of AO 0235+164 is about two times smaller 
than the critical accretion rate m CI ~ 0.02 in accretion disk evaporation model (Meyer et 
al. 2000; Liu et al. 2002) or much less than that of rh cr ~ 0.05 in strong ADAF principle 
(Narayan & Yi 1995; Abramowicz et al. 1995), implying that the central region of the 
accretion disk in AO 0235+164 is in ADAF (or RIAF) and the outer region of the accretion 
disk is a truncated standard thin disk. There is no consensus in literature about how to 
determine the transition radius R tT of inner thick ADAF and an outer thin standard disk. 
Three kinds of methods based on different physical processes have been proposed: disk 
evaporation (Meyer & Meyer-Hofmeister 1994; Liu et al. 1995b; Meyer et al. 2000; Liu 
et al. 2002; Roiahska & Czerny 2000), the absence of a strong ADAF solution (Narayan 
& Yi 1995; Abramowicz et al. 1995), and disk instability of inner radiation pressure 
dominated region with energy transfer outwards (Honma 1996; Kato & Nakamura 1998; Lu 
et al., 2004). The three methods generally give very different estimation of the transition 
radius. In the evaporation model, the transition radius at the critical accretion rate m cr is 
R tr ~ 355ro (Meyer et al. 2000; Liu et al. 2002), while in the strong ADAF scenario, the 
transition radius at the critical accretion rate m cr ~ 0.05 is about R tr — 560rc (Narayan & 
Yi 1995). The constraints from the existence of the broad line region in AGNs favored the 
relation of transition radius and accretion rate given by the strong ADAF principle (Czerny 



-19- 



et al. 2004). 

In the strong ADAF scenario with advection fraction of energy / = 0.9, the relation 
of the transition radius and accretion rate for rh < m cr is R tr ~ 0.074 (a/0.13) 4 m~ 2 r G 
(Narayan & Yi 1995), which gives a transition radius R tT — 1.5 x 10 3 r G . In the evaporation 
model (Meyer et al. 2000; Liu et al. 2002), 

R cv ~ 18.3m" a85 r G , (8) 

which gives a consistent transitional radius R tr = R cv — 1.2 x 10 3 r G . 

Rezzolla et al. (2003b) and Zanotti et al. (2003) showed that if it is globally perturbed, 
a torus or a thick disk of finite radial extent oscillates globally. The global oscillations lead 
to quasi-periodic accretion with frequencies of harmonic relationship 1:2:3:4. . . (Zanotti et 
al. 2003) 

fm = mfo, (9) 

where m — 1, 2, 3, . . .. The fundamental oscillation frequency fo is given with fo = c s /L 
where c s is sound speed and L = p/(dp/dr) is the radial length scale of pressure p. Note 
that / m in Zanotti et al. (2003) relates to the thick disk eigen-frequencies / n = [(2 + 
n)/2]/ for n = 0, 2, 4, . . . and m = (2 + n)/2 (Rezzolla et al. 2003b). Therefore, the 
trapped acoustic wave should have zero initial phase difference. The quasi-periodic p-mode 
oscillations of a small-size torus orbiting around a black hole have been suggested to explain 
the high-frequency QPOs (HFQPOs) in black hole X-ray binary systems (Rezzolla et al. 
2003a), which is observed to be in ratios of small integers (e.g. Abramowicz & Kluzniak 
2001). If the inner thick accretion flow of finite radial size in AO 0235+164 is perturbed by 
global perturbation or by strong local periodic agent, the accretion rate would vary quasi- 
periodically, triggering quasi-periodical plasma injection into relativistic jet due to jet-disk 
coupling and leading to the quasi-periodical variation of jet emission. As it describes an 
ADAF quite well except at two boundaries, the self-similar solution of ADAF is used in 
estimation of the fundamental oscillation frequency fo. From the self-similar solution for the 
ratio of specific heats 7 = 5/3 (Narayan & Yi 1994; Narayan et al. 2000), the pressure 

p oc pc 2 oc i?~ 5//2 and sound speed c s = 

a 2 < 1, where i>k is the Keplerian velocity at radius R. Then the fundamental frequency is 

fo = c s /L ~ (5/2)V2Q K 



, 1/2 

1 + ^) -1 



n l/2 

v K ~ (2/5) V 2 v K for 



- Lia^m . do) 



-3/2 

V.10 3 r G y V 1Q8M 

where Q K is the Keplerian angular frequency and r G = 2GM/c 2 = 2.95 x 10 13 cm (M/10 8 M Q ) 
is the Schwarzschild radius. The fundamental oscillation frequency of the thick ADAF of 
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AO 0235+164 is f ~ 0.182 yr^ 1 (R tT /1.2 x 10 3 r G ) _3/2 (M BH /4.72 x 1O 8 M ) _1 and the fun- 
damental period is P = 5.48yr (P tr /1.2 x 10 3 r G ) 3/2 (M BH /4.72 x 1O 8 M ). The oscillation 
model for torus or thick disk predicts a series of overtones in a sequence 1:2:3:. . ., independent 
of the disk radial extent and of the variation of sound speed with radius. The fundamental 
frequency and all the overtones are predicted to have nearly zero initial phase difference. 

With the redshift z = 0.94 of AO 0235+164, we have the expected lowest frequency 
/ obs = /o/(l + z) ~ 2.57 x lO-Mays- 1 ^/!^ x 10 3 r G )" 3/2 (M BH /4.72 x 1O 8 M ) _1 and 
period P obs ~ 10.7yr (P tr /1.2 x 10 3 r G ) 3/2 (M BH /4.72 x 1O 8 M ), respectively, which are con- 
sistent with the detected frequency v\ ~ (2.28 ± 1.10) x 10~ 4 days~ and period Pi ~ 
(12.0 ± 5.8) yrs (or 2P 2 = 10.9 yr at higher accuracy) at 8 GHz and other radio frequen- 
cies. As the observed lowest frequency, the harmonics and the difference of initial phases 
are consistent with the predictions of a thick ADAF oscillation model, we suggest that the 
multiple harmonic QPOs in the radio light curves of AO 0235+164 are due to the perturbed 
p-mode oscillations of an ADAF with a finite radial extent Pt r ~ 1.2 x 10 3 r G between a 
central supermassive black hole and an outer cool thin disk. 

There are two difficulties needed to be addressed with the model. The first one is the 
relative strength of QPOs in the power spectrum. In the theoretical power spectrum of 
oscillation accretion, the fundamental QPO is the strongest and the strength of harmonics 
decreases with the oscillation frequency (Rezzolla et al. 2003b; Zanotti et al. 2003). 
The Lomb power spectrum periodic analyses show that the rms of QPOs indeed decreases 
with frequency except the 6th harmonics with frequency uq = (1.511 ± 0.099) x 10~ 3 day _1 
or period Pq = 1.81 ± 0.12yr, which is exceptionally prominent and inconsistent with the 
prediction of the disk oscillation model. The other is on the nature of perturbation source. 
There is no obvious global or local perturbation source in a SMBH-thick disk system in AO 
0235+164. We will discuss the two difficulties in Sec. 6. 



6. Perturbation sources and supermassive black hole binary 

Although in Sec. 5.2 we suggest that the observed multiple harmonic QPOs in AO 
0235+164 would be due to the p-mode oscillations of a thick accretion disk (ADAF or RIAF), 
the nature of strong excitation mechanisms, either local (Rubio-Herra & Lee 2005a,b) or 
global (Zanotti et al. 2003; Rezzolla et al. 2003a), is still an open question, in particularly 
for a system of a thick disk with large radial extent around a SMBH. Here we suggest that 
a SMBHB embedding in the thick accretion disk would serve as such a local, strong and 
periodic perturbation source as required by the p-mode oscillations of a thick accretion disk 
investigated by Rubio-Herra & Lee (2005b). A SMBHB may form in galaxy minor merger 
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and is introduced in a helical jet model for AO 0235+164 (Ostorero et al. 2004) and later in 
explaining the rapid changes in the VLBI (Very Long Baseline Interferometry) jet position 
angles of the radio source (Frey et al. 2006). A SMBHB embedding in an ADAF or RIAF 
cannot open a gap in an accretion disk, because to open a gap the disk opening angle H/R 
should be H/R < 0.3 and the mass ratio q of the secondary and the primary black holes 
should be q < 1 and q > q min = (8br/ '8) a (H/R) 2 > 3 for an ADAF with H/R ~ 1 and 
0.1 < a < 0.3 (Lin & Papaloizou 1986). The interaction with an RIAF is hydro-dynamically 
unimportant to the evolution of a SMBHB (Liu 2004), but the gravitational attraction of 
the secondary to the plasma in the thick disk changes periodically and triggers a series 
of density waves which propagate inwards and outwards in the thick disk with a limited 
radial extent. When they arrive at the inner or the outer boundaries, the acoustic waves 
are reflected backwards with a very large fraction of the incoming energy, if the acoustic 
frequency is about the fundamental frequency or harmonics of the cavity. Other acoustic 
waves are almost completely be absorbed by the boundaries. This kind of p-mode oscillations 
in a tours or thick accretion disk were first discussed for global perturbation by Rezzolla et 
al. (2003b) and later for local strong periodic agent by Rubio-Herra & Lee (2005a). A 
SMBHB thus perfectly serves as a local strong perturbation agent to a SMBH-thick disk 
system. 

For a SMBHB with separation a < R tr , the orbital period Pb should be 

1/2 



Pb = 2tt 



a 3 



GM BH (1 + q) 



Rt 



3/2 



In the observer's frame, the period is observed to be P G bs = (1 + z)Pb < 80 yr, which may 
or may not be within the detection window of the three databases which we used. Rubio- 
Herra & Lee (2005b) showed that in addition to the strong QPO with very high quality 
factor Q at the perturbation period Pb, the oscillations of a perturbed thick disk should be 
detected significantly at Pb/2 and possibly also at 2P b . In the helical-jet model (Ostorero 
et al. 2004), the SMBHB orbital period is identified with P obs = P 2 ~ 5 - 6yrs. If the 
helical jet model is correct, the SMBHB has a separation a ~ 166rc, which is much smaller 
than the transitional radius R tr and implies an embedding SMBHB in ADAF (or RIAF). In 
this scenario, SMBHB should excite thick disk oscillations not only at period P 2 but also at 
periods P 4 and probably P 1 . One difficulty with the explanation is that one has to explain 
the exceptional prominence of the 6th QPO with very high quality factor Q, which should 
be the weakest among six QPOs and may not be detectable in the disk oscillation model. 
Therefore, a more reasonable suggestion would be that the 6th QPO of period P 6 = 1.81 yr 
corresponds to the SMBHB orbital motion. In this scenario, a QPO at period Pe/2 = 0.91 yrs 
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should also be significant, which might have already been detected with P = 1.03 ± 0.03 yrs 
at the Lomb power spectrum at 8 GHz in Fig. 2. However, a QPO with a period about one 
year may not be identified correctly due to the one-year astronomical cycle. 



For a SMBHB with an orbital period P obs = 1.81 yr or 5.46 yr, the binary separation is 

-fobs 



a = r G 



2/3 

(l+g) 1/3 



1 + z)2^2nr G _ 
~ 166r G (1 + q) 1/3 if P obs = 5.46 yr, or 

~ 79r G (l + g) 1/3 if P obs = 1.81 yr. (12) 

Because the interaction with an ADAF or RIAF is negligible and the evolution of SMBHB 
is dominated by the gravitational wave radiation (Liu 2004), the lifetime of SMBHB with a 
circular orbit in AO 0235+164 due to the gravitational wave radiation is (Peters & Mathews 
1963) 

a 5 / a \ 4 i . . 
*gw = — =o( — J + 



a\ 8 \r G 

~ 7.0 x 10V 1 (1 + # /3 if Pobs = 5.46 yr, or 

~ 3.6x 10 5 yr (^^'(l + g) 173 if P obs = 1.81 yr. (13) 

Both identifications of the SMBHB orbital motion suggest a short-lived binary. However, it 
is also possible that the separation of SMBHB is a ~ 10 3 r G and the orbital period P b ~ 70 yr 
is out of the observational window. A SMBHB with such a large separation is long-lived with 
t gw ~ 9.2 x 10 9 yr (q^-) (1 + qf^ ■ The observations of the changes of VLBI jet positional 
angle may help to resolve the puzzles. 



7. Discussions and conclusions 

We have investigated the radio variabilities of AO 0235+164 and analyzed the periods 
of radio light curves at 4.8, 8, 14.5, 22, and 37 GHz, basing on the databases of UMRAO, 
NRAO and Metsahovi Observatory. We first study the consistence of the UMRAO 's and 
NRAO's data at 8 GHz for AO 0235+164 and then construct a combined radio light curve 
with the observational data from the two databases. The periodic analyses with three kinds 
of classical periodic analysis methods show with unprecedentedly high signal-to-noise ratio 
that the variations of the combined radio light curves are composed of six periodic variations 
plus random fluctuations. In the power spectrum, the peaks are Gaussian and identified with 
QPOs of periods Pi = 12.02 yr, P 2 = 5.45 yr, P 3 = 3.62 yr, P 4 = 2.83 yr, P 5 = 2.15 yr, and 
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P 6 = 1.81 yr. The coherent quality factor Q of the six QPOs are limited by the monitoring 
program time. We discover a harmonic relationship for the six QPOs, based on the ratio of 
QPO frequencies in a sequence of small integers 1:2:3:4:5:6 and on the zero or n difference 
of initial phases of QPOs. The harmonic relationship suggests that the six QPOs should 
have the same physical origin. Among the six harmonic QPOs, the second one with period 
P 2 ~ 5.46 ± 0.47 yrs is the strongest and is identified with the fundamental period. The five 
strong harmonics reported in some literature are confirmed here with an unprecedentedly 
high signal-to-noise and accuracy, while the weakest QPO of period P 5 = 2. 15 yrs is reported 
first time in this paper. 

We investigate the dependence of harmonic relationship, coherent quality factor, the 
relative rms, the differences of the initial phase of QPOs on radio frequency. The results 
show that all the QPO properties are independent of radio frequencies. The six harmonic 
QPOs give a good fit to the outburst structures of the radio light curves at 4.8, 8, 14.5, 22 
and 37 GHz, and predict major outbursts and outburst structures significantly different from 
that given with a single period of P ~ 5.46 yr (Raiteri et al. 2001). Both the single-period 
and the multiple harmonic QPO (for the low frequencies) scenarios predict a major outburst 
in the early half of 2004, while the combinations of six harmonic QPOs suggest only minor 
fluctuations during that period but a major burst in the late half of 2004 or the early half 
2005 at 37 GHz, as it is shown in Fig. 7. In our periodic analyses of radio light curves at 22 
GHz and 37 GHz, we have used the data after 2001 and the harmonic QPO scenario predict 
a major outburst in the early 2005. At lower frequencies we cannot get the observational 
data in recent years and make the analyses basing on the data before 2001. Therefore, the 
predictions are slightly different for different wave-bands and the predicted major bursts 
may be in the early 2005. However, the observations of a WEBT (the Whole Earth Blazar 
Telescope) campaign started from 2003 show that AO 0235+164 is quiet in recent years and 
the predicted major outbursts in the spring of 2004 did not come (Raiteri et al. 2005). 
Is the predicted major outburst just delayed significantly or is the activity of the object 
interrupted due to some reasons? When will the object re-activate with a different active 
cycle? To answer the questions, more observations for another several years are needed. 

Either a significant delay of outbursts or the interruption of activity for some time 
would exclude the single period scenario and the simple version of a helical jet model but 
may be consistent with the SMBHB-thick disk oscillation model. In a SMBHB-thick disk 
oscillation scenario, the multiple QPOs in the radio light curves of AO 0235+164 are due 
to quasi-periodic injection of plasma from an oscillation thick disk of a finite radial extent 
into relativistic jet. A torus or a thick accretion disk of finite radial extent oscillates and 
accretes matter quasi-periodically with small integer ratio of frequencies with 1:2:3:. . ., if it 
is perturbed by a global perturbation or a local strong periodic agent (Zanotti et al. 2003; 
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Rezzolla et al. 2003a; Rubio-Herra & Lee 2005b). The fundamental oscillation frequency 
of a thick disk is determined by the radial extent and sound speed. We estimated the 
central black hole mass Mbh = (4.72 ± 2.04) x 1O 8 M and the dimensionless accretion rate 
fn ~ 0.007 of AO 0235+164, based on the observations of emission line properties. The 
estimated dimensionless accretion rate with the current knowledge of accretion around black 
hole suggests that the inner region of accretion disk in AO 0235+164 is a geometrically thick 
disk (i.e. ADAF or RIAF) with radial extent about 1.2 x 10 3 r G and the outer region is a 
truncated cool thin disk. With the self-similar solution of ADAF, we analytically estimate the 
fundamental oscillation frequency / Q b s — 0.0938 yr _1 or period P obs = 10.7 yr in the observer's 
frame, which together with the predicted overtones are consistent with the observed multiple 
QPOs with lowest frequency P\ ~ 12.02 ±5.79yr (or 2P 2 = 10.9 yrs at higher accuracy) and 
harmonics. As there is no natural global or strong local perturbation source in standard 
SMBHB-accretion disk system to excite the global oscillations of a thick disk of a large 
radial extent Rd ~ 10 3 rc, we suggest a SMBHB rotating within the thick disk to excite the 
p-mode oscillations of RIAF in AO 0235+164. The rotation of the secondary triggers acoustic 
waves in an ADAF (or RIAF), which propagate inwards and outwards in the thick disk and 
are reflected backwards by the inner and outer boundaries. The p-mode oscillations with 
frequencies of about the fundamental oscillation frequency and the overtones are trapped 
by the cavity. The orbital motion of a SMBHB could lead to the formation of a helical 
morphology of relativistic jet and to the rapid changes of VLBI jet position angles. Therefore, 
in the SMBHB-disk oscillation model, the relativistic jet can be also helical. The orbital 
motion may correspond to the 6th QPO with P 6 = 1.81 yrs, the 2nd QPO with P 2 = 5.46 yrs, 
or have not been observed. Long term VLBI monitoring of rapid changes of jet positional 
angles would help to answer the question. One of the differences with the simple helical jet 
model is that in the SMBH-thick disk oscillation scenario the minor burst events in the radio 
light curves is no longer random but together with the major outbursts are the combinations 
of harmonic six QPOs which are the p-mode oscillations of the RIAF and the binary orbital 
motion. This is a first time report of multiple harmonic QPOs and the possible detection of 
perturbed disk p-mode oscillations in AGNs. 

The second difference of predictions of the helical jet model and the SMBHB-thick disk 
oscillation scenario is that periodic outbursts in the former model are persistent and regular 
while in the later scenario are temporary with moderate coherence quality factor Q and 
depend on the dimensionless accretion rate. In the SMBHB-thick oscillation scenario, reso- 
nant harmonic QPOs would disappear from the radio light curves or the periodic activities 
of the radio source are interrupted, when the accretion rate increases only about two times 
aand becomes larger than the critical accretion rate m cr and the inner region of accretion 
disk becomes a geometrically thin cool disk. Therefore, the quiescence and the absence of 
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outbursts from the radio light curves of AO 0235+164 may be because of the transition of 
accretion mode from a RIAF to a standard thin cool accretion disk due to the moderate 
increase of accretion rate. The X-ray spectral observations with the Chandra satellite in 
August 2000 and with the XMM-Newton satellite in August 2004 show tentative detections 
of the red-shifted Fe Kq fluorescent emission line in AO 0235+164 (Raiteri et al. 2006), 
which suggests that the accretion disk at the immediate vicinity of the SMBH is very likely a 
standard thin accretion disk (Fabian et al. 2000). Therefore, the absence of outbursts from 
the radio light curves and the interruption of the harmonic QPOs can be understood in the 
SMBHB-thick disk oscillation scenario as the transition of accretion mode due the increase 
of accretion rate. In a jet-disk coupling scenario, an object becomes radio quiet when the 
accretion disk becomes a standard thin disk (Fender & Belloni 2004). 

We notice that in the helical jet model, the separation of SMBHB is much smaller than 
the transitional radius R tr of inner thick RIAF and outer thin disk. Therefore, the SMBHB 
in the helical jet model for AO 0235+164 should also excite the p-mode oscillations of an 
RIAF, leading to quasi-periodic accretion and quasi-periodic injection of plasma into the 
relativistic jet. Therefore, multiple QPOs with harmonic relationship are also expected by 
the helical jet model for AO 0235+164, if we consider the RIAF p-mode oscillations. 

The embedding secondary black hole spherically accretes matter from the RIAF and 
radiates energy in X-ray. The accretion radius of the secondary black hole is r acc gjr^p , 
where Av is the differential velocity of the secondary black hole and the plasma in the thick 
disk. For a minor merger with g C 1, we obtain r acc <C R and (Av) 2 <C c 2 . From the 
self-similar solution c 2 = (2/5) V£ for 7 = 5/3 (Narayan & Yi 1994; Narayan et al. 2000), 
we have r acc (5/2)g(l + q)R. The accretion rate of the secondary black hole is 

M s ~ 4vrr 2 cc c s p, (14) 

where the mass density p is given with the self-similar solution and p = 1 2?/| 7r ^Edd^o: - 1 R~ 2 Vk 1 • 
From Eq. (14), we have M s ~ 2 ^y^ g 2 (l + q) 2 M^^ma 1 . Defining the Eddington accretion 
rate for the secondary black hole M| dd = L| dd /0.1c 2 with L| dd = 1.26 xlO 38 (m/M ) ergs s _1 , 
we have the dimensionless accretion rate of the secondary black hole, m s = M s /M^ dd ~ 
6.6(1 + q) 2 (q/a)m. For the typical parameters q = 0.01 and a = 0.13, the relative accretion 
rate of the secondary is m s ~ 0.51m ~ 5.1 x 10 -3 , while m s ~ 5.1 x 10~ 2 for q = 0.1. The 
X-ray emission from the accretion of the secondary black hole is small and may not be ob- 
servable in AO 0235+164, because of the strong X-ray radiations both from the relativistic 
jet and from accretion disk around the primary. A different 7 value does not change the 
conclusion. 



Our work implies that if a thick accretion disk with finite radial extent is perturbed 
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globally or locally by, e.g. SMBHB, the p-mode oscillations can be observed as multiple 
harmonic QPOs, when the relativistic jet is along the line of sight. Periods have been reported 
and SMBHBs have been suggested in literature for many AGNs. So, multiple QPOs with 
a harmonic relationship of frequencies are expected by the SMBHB-disk oscillation model 
in moderately low luminosity AGNs with SMBHB at center, which may have accretion 
configuration like the one in the BL Lac object AO 0235+164. We will report our periodic 
analysis results on a large sample of AGNs in a coming paper (Zhao, Liu, & Wu 2006). 
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Fig. 3. — The power spectrum of the radio light curves at 8.0 GHz at UMRAO (left) and 
at 8.2 GHz at NRAO (right). The identified QPOs are indicated. The dashed and the 
dash-dotted lines are the la and 3a significance levels, respectively. 
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Fig. 4. — Quality factor Q vs repetition number N of QPOs in the combined radio light 
curves. The correlation of Q and N implies that the Q-value is limited by the monitoring 
time of the object and much longer observations are needed to give stringent constraint on 
the intrinsic coherence of QPOs. 
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Fig. 5. — Normalized Jurkevich periodic analysis results of the combined radio light curve 
at 8 GHz. The minima corresponding to periods 1.78 yrs, 3.46 yrs, 5.34 yrs, 7.90 yrs and 
11.60 yrs are prominent. A period is strong if / > 0.25 with / = (1 — V^j/V^ 
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Fig. 6. — The periodic analysis results with ZDCF method. The results are consistent with 
that of Jurkvich method and the two weakest QPOs at Lomb's periodogram are missing. 
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Fig. 7. — The radio light curves of AO 0235+164. The observation frequencies are indicated 
at the upper right corner in the figures. The solid line is the result fitted with six QPOs 
and the dotted line is the result fitted with the fundamental QPO of period P s . The results 
clearly showed that in addition to the fundamental QPO the other five QPO harmonics 
significantly contributed to the variations of the radio light curves. 
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Fig. 8. — Lomb periodograms of the radio light curves at 4.8 GHz (uper left), 14.5 GHz 
(uper right), 22 GHz (lower left), and 37 GHz (lower right). Fitted centroid periods of QPOs 
with Gaussian function are indicated. The dashed and dash-dotted lines are the la and 3a 
significance levels, respectively. The radio data at 4.8 GHz and 14.5 are from UMRAO and 
at 22 GHz and 37 GHz from Metsahovi Observatory. 
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Fig. 9. — Dependence of QPO frequencies on radio frequencies. The horizontal dashed lines 
with small integer number N = 1,2, 3, 4, 5, and 6 indicates the expected ratio of frequencies 
of QPOs and the basic QPO with lowest frequency (calculated basing the fundamental QPO 
frequency u s at 8 GHz). The centroid frequencies and the harmonic relationship of QPOs 
are independent of the observational radio frequencies. 
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Fig. 10. — The amplitude rms (in percentage) of QPO harmonics as a function of the radio 
frequencies {upper panel) and the amplitude rms relative to that of the fundamental QPO 
(the 2nd harmonics) vs radio frequencies (lower panel). The solid line is the fitted results 
with least square method. 
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Table 1: Hormanic QPOs in the radio light curves at 4.8 GHz, 8 GHz, 14.5 GHz, 22 GHz 
and 37 GHz analyzed with the Lomb power spectrum method. 
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1.81±0.12 


12.98 


0.332 


14.5 


3.00 


6 


0.07 


4.8 GHz 


0.2568±0.0797 


10.86±0.66 


10.67±3.31 


2.74 


0.132 


2.3 


0.54 


1 


-0.35 




0.4716±0.0700 


9.02±0.79 


5.81±0.86 


5.72 


0.031 


4.2 


1.00 


2 


0.00 


5.62 


0.7254±0.0771 


8.20±0.87 


3.78±0.40 


7.99 


0.085 


6.5 


1.54 


3 


-0.19 


±0.44 


0.9295±0.0780 


7.08±1.01 


2.95±0.25 


10.12 


0.019 


8.5 


1.97 


4 


1.19 


















(5) 


(-0.41) 




1.522±0.066 


8.27±0.86 


1.80±0.08 


19.53 


0.006 


13.6 


3.23 


6 


0.36 


14.5 GHz 


0.2019±0.0558 


6.19±0.57 


13.57±3.75 


3.07 


0.014 


1.9 


0.42 


1 


-0.34 




0.4826±0.0532 


8.10±0.44 


5.68±0.63 


7.70 


0.251 


4.4 


1.00 


2 


0.00 


5.53 


0.7299±0.0516 


6.21±0.57 


3.75±0.27 


12.01 


0.033 


6.7 


1.51 


3 


-0.11 


±0.33 


0.9776±0.0568 


7.58±0.47 


2.80±0.16 


14.62 


0.041 


9.0 


2.03 


4 


1.22 




1.297±0.060 


4.95±0.72 


2.11±0.10 


18.36 


0.002 


11.9 


2.69 


5 


-0.40 




1.472±0.067 


6.72±0.53 


1.86±0.08 


18.77 


0.004 


13.5 


3.05 


6 


0.44 


22 GHz 


0.3138±0.0625 


7.73±0.23 


8.73±1.74 


4.26 


0.087 


2.6 


0.66 


1 


-0.12 




0.4755±0.0924 


8.96±0.20 


5.76±1.12 


4.37 


0.087 


3.9 


1.00 


2 


0.00 


5.55 


0.6890±0.0701 


7.47±0.23 


3.98±0.40 


8.35 


0.025 


5.6 


1.45 


3 


-0.34 


±0.47 


1.027±0.079 


6.06±0.29 


2.67±0.21 


11.04 


0.015 


8.4 


2.16 


4 


0.94 


















(5) 


(-0.29) 




1.516±0.078 


6.56±0.27 


1.81±0.09 


16.59 


0.015 


12.4 


3.19 


6 


-0.02 


37 GHz 


0.1214±0.0753 


5.99±0.27 


22.57±14.00 


1.37 


0.039 


1.03 


0.30 


1 


-0.55 




0.4065±0.0948 


7.07±0.23 


6.74±1.57 


3.64 


0.040 


3.5 


1.00 


2 


0.00 


5.77 


0.7376±0.0644 


3.95±0.41 


3.71±0.32 


9.73 


0.040 


6.3 


1.81 


3 


0.19 


±0.64 


0.9879±0.1305 


4.75±0.41 


2.78±0.37 


10.74 


0.073 


8.4 


2.43 


4 


-0.42 


















(5) 


( 0.36) 




1.432±0.113 


5.93±0.27 


1.91±0.15 


9.13 


0.073 


12.2 


3.52 


6 


-0.77 



